rm(list = ls())
library(deSolve)

n = 0;
rho = 0.02;
delta = 0.05;
theta = 3;
alpha = 0.36;
A = 5;

# 线性化的微分方程才行，下面是非线性化的，所以不一定在鞍形稳定路径上，
# 存在发散行为
ramsey <- function(t, y, parms){
  with(as.list(y),{
    dk <- A*k^alpha - cons -(n + delta)*k
    dcons <- (cons/theta) * (A*alpha*k^(alpha-1)-delta-rho)
    list(c(dk,dcons))
  })
}

kss <- ((rho + delta)/(A*alpha))^(1/(alpha-1))
css <- A*kss^alpha-(n+delta)*kss

tt <- seq(0,50,1)
out <- ode(y = c(k = kss-20,cons = css-1), times = tt, func = ramsey, parms = NULL, method = rkMethod('rk45ck'))
picdata <- as.data.frame(out)
# ggplot(picdata, aes(x = k, y = cons)) + geom_line() + theme_bw()
plot(out)

